{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse.linalg import spsolve\n",
    "import scipy.sparse as sps\n",
    "import sympy\n",
    "from math import pi\n",
    "from time import time\n",
    "\n",
    "import pyamg\n",
    "\n",
    "import setup_grids\n",
    "\n",
    "from porepy.numerics.fv import mpsa, fvutils\n",
    "\n",
    "from porepy.params.tensor import FourthOrderTensor as stiffness\n",
    "from porepy.params import bc\n",
    "from porepy.viz.plot_grid import plot_grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.011010878448828448, 0.0030197500225270933, 0.0007837446628904635, 0.00019841406340054548] displacement error\n",
      "[0.036391713833312, 0.010752498606472382, 0.003027389815393291, 0.0008407990967767796] stress error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[0.013141424612606898, 0.003959554975380245, 0.001031705289450405, 0.00024384620727654387] displacement error\n",
      "[0.03618728579049561, 0.015764291256799894, 0.005387936243555718, 0.0020697006664011576] stress error\n",
      "Triangular errors - no perturbations\n",
      "[0.007499534100939774, 0.002035510853495113, 0.0005339430679082081, 0.0001366067446813222] displacement error\n",
      "[0.048716801365141, 0.015930725831506523, 0.005076956913423936, 0.0016467565321741852] stress error\n",
      "Triangular errors - perturbation 0.5\n",
      "[0.009376925921720678, 0.0026191763163007998, 0.0006930802706380639, 0.00017435555727363358] displacement error\n",
      "[0.04808796652224495, 0.019133502522555433, 0.007100574505058305, 0.00277664388741696] stress error\n"
     ]
    }
   ],
   "source": [
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "lmbda = 1\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 4\n",
    "grid_type = 'cart'\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "ux = sympy.sin(x) * sympy.cos(y)\n",
    "uy = sympy.sin(x) * x**2\n",
    "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "divu = dux_x + duy_y\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n",
    "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "            \n",
    "        # Permeability tensor\n",
    "        mu_c = mu * np.ones(g.num_cells)\n",
    "        lmbda_c = lmbda * np.ones(g.num_cells)\n",
    "        k = stiffness(2, mu_c, lmbda_c)\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.get_all_boundary_faces()\n",
    "        bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        stress, bound_stress = mpsa.mpsa(g, k, bound_cond)\n",
    "        div = fvutils.vector_divergence(g)\n",
    "        a = div * stress\n",
    "        \n",
    "        # Boundary conditions\n",
    "        xf = g.face_centers\n",
    "        u_bound = np.zeros((g.dim, g.num_faces))\n",
    "        u_bound[0, bound_faces] = ux_f(xf[0, bound_faces], xf[1, bound_faces])\n",
    "        u_bound[1, bound_faces] = uy_f(xf[0, bound_faces], xf[1, bound_faces])\n",
    "        \n",
    "        \n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n",
    "        b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        ux_num = u_num[::2]\n",
    "        uy_num = u_num[1::2]\n",
    "        \n",
    "        stress_x_num = stress_num[::2]\n",
    "        stress_y_num = stress_num[1::2]\n",
    "        # Exact solution\n",
    "        ux_ex = ux_f(xc[0], xc[1])\n",
    "        uy_ex = uy_f(xc[0], xc[1])\n",
    "        u_ex = np.vstack((ux_ex, uy_ex))\n",
    "        u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "        \n",
    "        sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n",
    "        sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n",
    "        \n",
    "        stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n",
    "        stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n",
    "        stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                 stress_y_num - stress_y_ex))\n",
    "        stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n",
    "        \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * u_diff**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(f_cart_nopert,'stress error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'displacement error')\n",
    "print(f_cart_pert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(f_triang_nopert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('triangular', 0.5)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'displacement error')\n",
    "print(f_triang_pert,'stress error')\n",
    "\n",
    "# EK: These values were hard-coded 12.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] - 0.00078374466288952144) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.0030273898153929924) < 1e-10\n",
    "assert np.abs(u_cart_pert[2] - 0.0010317052894504846) < 1e-10\n",
    "assert np.abs(f_cart_pert[2] - 0.0053879362435559708) < 1e-10\n",
    "\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018. The values changed because\n",
    "# the default value for eta for triangular grids was changed to 1/3\n",
    "assert np.abs(u_triang_nopert[2] - 0.0005339430679082081) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.005076956913423936) < 1e-10\n",
    "assert np.abs(u_triang_pert[2] - 0.0006930802706380639) < 1e-10\n",
    "assert np.abs(f_triang_pert[2] - 0.007100574505058305) < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neumann boundary conditions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.19341322457532645, 0.05398606173752483, 0.019785099166683015] displacement error\n",
      "[3.58265112 2.72862225] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.03080448569971605, 0.022727907422006446, 0.017633923832997576] stress error\n",
      "[1.35535952 1.28887408] convergence rate f^n-1 / f^n\n",
      "\n",
      "Triangular errors - no perturbations\n",
      "[0.15803686415889875, 0.04310823056510322, 0.015475821210292553] displacement error\n",
      "[3.6660485  2.78552136] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.04467366719714846, 0.02540025474312989, 0.016373893423034032] stress error\n",
      "[1.75878816 1.55126543] convergence rate f^n-1 / f^n\n"
     ]
    }
   ],
   "source": [
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "lmbda = 1\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 3\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "ux = (0.5-y)**2 * sympy.sin(y) * sympy.cos(y)\n",
    "uy = y * (y-1) * sympy.cosh(x)\n",
    "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "divu = dux_x + duy_y\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n",
    "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "            \n",
    "        xf = g.face_centers\n",
    "        sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n",
    "        sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n",
    "        \n",
    "        stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n",
    "        stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n",
    "\n",
    "        # Permeability tensor\n",
    "        mu_c = mu * np.ones(g.num_cells)\n",
    "        lmbda_c = lmbda * np.ones(g.num_cells)\n",
    "        k = stiffness(2, mu_c, lmbda_c)\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        n = g.nodes.max()\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        top = np.ravel(np.argwhere(g.face_centers[1, :] > n - 1e-10))\n",
    "        bot = np.ravel(np.argwhere(g.face_centers[1, :] < 1e-10))\n",
    "        left = np.ravel(np.argwhere(g.face_centers[0, :] < 1e-10))\n",
    "        right = np.ravel(np.argwhere(g.face_centers[0, :] > n - 1e-10))\n",
    "\n",
    "        dir_faces = np.hstack((bot, top))\n",
    "        neu_faces = np.hstack((left, right)) \n",
    "        bound_cond = bc.BoundaryCondition(g, dir_faces, ['dir'] * dir_faces.size)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        stress, bound_stress = mpsa.mpsa(g, k, bound_cond)\n",
    "        div = fvutils.vector_divergence(g)\n",
    "        a = div * stress\n",
    "        xf = g.face_centers\n",
    "        \n",
    "        nfi, _, sgn = sps.find(g.cell_faces[neu_faces,:])\n",
    "        \n",
    "        u_bound = np.zeros((g.dim, g.num_faces))\n",
    "        u_bound[0, dir_faces] = ux_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        u_bound[1, dir_faces] = uy_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        \n",
    "        # When setting the Neumann boundary condition we need to flip the \n",
    "        # sign of the faces with a normal pointing inwards as it is assumed\n",
    "        # that the Neumann condition is the force from the boundary on to face\n",
    "        u_bound[0, neu_faces[nfi]] = stress_x_ex[neu_faces[nfi]] * (sgn)\n",
    "        u_bound[1, neu_faces[nfi]] = stress_y_ex[neu_faces[nfi]] * (sgn)\n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n",
    "        b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        ux_num = u_num[::2]\n",
    "        uy_num = u_num[1::2]\n",
    "        \n",
    "        stress_x_num = stress_num[::2]\n",
    "        stress_y_num = stress_num[1::2]\n",
    "        # Exact solution\n",
    "        ux_ex = ux_f(xc[0], xc[1])\n",
    "        uy_ex = uy_f(xc[0], xc[1])\n",
    "        u_ex = np.vstack((ux_ex, uy_ex))\n",
    "        u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "        \n",
    "        stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                 stress_y_num - stress_y_ex))\n",
    "        stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n",
    "        \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * u_diff**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2)))\n",
    "                \n",
    "    return u_err, flux_err\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "u_cart_rate = np.array(u_cart_nopert)\n",
    "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
    "f_cart_rate = np.array(f_cart_nopert)\n",
    "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_cart_nopert,'stress error')\n",
    "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "u_triang_rate = np.array(u_triang_nopert)\n",
    "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
    "f_triang_rate = np.array(f_triang_nopert)\n",
    "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
    "\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_triang_nopert,'stress error')\n",
    "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
    "\n",
    "# # EK: These values were hard-coded 12.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] - 0.01978509916668196) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.017633923832997497) < 1e-10\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2] - 0.015475821210292553) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.016373893423034032) < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robin Condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.13785825178140515, 0.033672908859893956, 0.0082926160732301] displacement error\n",
      "[4.09404047 4.06058939] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.04932864932955223, 0.01582468790884128, 0.004760957079948613] stress error\n",
      "[3.11719571 3.32384595] convergence rate f^n-1 / f^n\n",
      "\n",
      "Triangular errors - no perturbations\n",
      "[0.11226477176063937, 0.025455746327683503, 0.006007054619672307] displacement error\n",
      "[4.41019369 4.2376419 ] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.044251086212033885, 0.016626084791879308, 0.0058474670379793795] stress error\n",
      "[2.6615458  2.84329688] convergence rate f^n-1 / f^n\n"
     ]
    }
   ],
   "source": [
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "lmbda = 1\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 3\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "ux = (0.5-y)**2 * sympy.sin(y) * sympy.cos(y)\n",
    "uy = y * (y-1) * sympy.cosh(x)\n",
    "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "divu = dux_x + duy_y\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n",
    "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "            \n",
    "        xf = g.face_centers\n",
    "        sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n",
    "        sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n",
    "        \n",
    "        stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n",
    "        stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n",
    "\n",
    "        # Permeability tensor\n",
    "        mu_c = mu * np.ones(g.num_cells)\n",
    "        lmbda_c = lmbda * np.ones(g.num_cells)\n",
    "        k = stiffness(2, mu_c, lmbda_c)\n",
    "        robin_weight = 1.0\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        n = g.nodes.max()\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        top = np.ravel(np.argwhere(g.face_centers[1, :] > n - 1e-10))\n",
    "        bot = np.ravel(np.argwhere(g.face_centers[1, :] < 1e-10))\n",
    "        left = np.ravel(np.argwhere(g.face_centers[0, :] < 1e-10))\n",
    "        right = np.ravel(np.argwhere(g.face_centers[0, :] > n - 1e-10))\n",
    "\n",
    "        dir_faces = np.hstack((top, left))\n",
    "        rob_faces = np.hstack((right, bot))\n",
    "        bc_faces = np.hstack([dir_faces, rob_faces])\n",
    "        bc_names = ['dir']*len(dir_faces) + ['rob'] * len(rob_faces)\n",
    "        bound_cond = bc.BoundaryCondition(g, bc_faces, bc_names)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        stress, bound_stress = mpsa.mpsa(g, k, bound_cond,robin_weight=robin_weight)\n",
    "        div = fvutils.vector_divergence(g)\n",
    "        a = div * stress\n",
    "        xf = g.face_centers\n",
    "        \n",
    "        rfi, _, sgn = sps.find(g.cell_faces[rob_faces,:])\n",
    "        \n",
    "        u_bound = np.zeros((g.dim, g.num_faces))\n",
    "        u_bound[0, dir_faces] = ux_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        u_bound[1, dir_faces] = uy_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        \n",
    "        # When setting the Neumann boundary condition we need to flip the \n",
    "        # sign of the faces with a normal pointing inwards as it is assumed\n",
    "        # that the Neumann condition is the force from the boundary on to face\n",
    "        u_bound[0, rob_faces[rfi]] = stress_x_ex[rob_faces[rfi]] * (sgn) + \\\n",
    "            robin_weight * ux_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n",
    "\n",
    "        u_bound[1, rob_faces[rfi]] = stress_y_ex[rob_faces[rfi]] * (sgn) +\\\n",
    "            robin_weight * uy_f(xf[0, rob_faces[rfi]], xf[1, rob_faces[rfi]]) * g.face_areas[rob_faces[rfi]]\n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n",
    "        b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        ux_num = u_num[::2]\n",
    "        uy_num = u_num[1::2]\n",
    "        \n",
    "        stress_x_num = stress_num[::2]\n",
    "        stress_y_num = stress_num[1::2]\n",
    "        # Exact solution\n",
    "        ux_ex = ux_f(xc[0], xc[1])\n",
    "        uy_ex = uy_f(xc[0], xc[1])\n",
    "        u_ex = np.vstack((ux_ex, uy_ex))\n",
    "        u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "        \n",
    "        stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                 stress_y_num - stress_y_ex))\n",
    "        stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n",
    "        \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * u_diff**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2)))\n",
    "                \n",
    "        \n",
    "    return u_err, flux_err\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "u_cart_rate = np.array(u_cart_nopert)\n",
    "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
    "f_cart_rate = np.array(f_cart_nopert)\n",
    "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_cart_nopert,'stress error')\n",
    "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "u_triang_rate = np.array(u_triang_nopert)\n",
    "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
    "f_triang_rate = np.array(f_triang_nopert)\n",
    "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
    "\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_triang_nopert,'stress error')\n",
    "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
    "\n",
    "# RB: These values were hard-coded 11.09.2018.\n",
    "assert np.abs(u_cart_nopert[2] - 0.0082926160732301) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.004760957079948613) < 1e-10\n",
    "assert np.abs(u_triang_nopert[2] - 0.006007054619672307) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.0058474670379793795) < 1e-10\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that the reduction rates are smallar than for the pure Dirichlet case - the asymptotic rates are 1 for displacement and 0.5 for stress. The reason for this is not understood - there are no obvious bugs in the code. Todo.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[[0.37388856 0.4014268  0.38285806]\n",
      " [0.08186136 0.09736873 0.08866018]\n",
      " [0.01980355 0.02491816 0.02207893]] displacement error\n",
      "[[0.16987677 0.16287318 0.17095043]\n",
      " [0.06445937 0.05898811 0.06587345]\n",
      " [0.02181871 0.01941263 0.02240798]] stress error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[[0.38707647 0.51439666 0.52211052]\n",
      " [0.10514942 0.11284306 0.10794662]\n",
      " [0.02294892 0.02818847 0.02665215]] displacement error\n",
      "[[0.22785554 0.22643673 0.23802875]\n",
      " [0.07617105 0.06993615 0.07181565]\n",
      " [0.02941291 0.02754268 0.02920611]] stress error\n",
      "Triangular errors - no perturbations\n",
      "[[0.28239685 0.28600247 0.27724259]\n",
      " [0.06450559 0.06808808 0.0647828 ]\n",
      " [0.01574973 0.01710995 0.01599735]] displacement error\n",
      "[[0.14821803 0.14406458 0.14628784]\n",
      " [0.05653494 0.05502727 0.05612567]\n",
      " [0.01890709 0.01837107 0.01901593]] stress error\n",
      "Triangular errors - perturbation 0.5\n",
      "[[0.34369963 0.45755484 0.47312797]\n",
      " [0.09035705 0.09202687 0.08932667]\n",
      " [0.01911508 0.02070007 0.02009385]] displacement error\n",
      "[[0.18615258 0.18026657 0.17874808]\n",
      " [0.06705379 0.06574299 0.06617049]\n",
      " [0.02295011 0.02251892 0.02313161]] stress error\n"
     ]
    }
   ],
   "source": [
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "pert = 0.5\n",
    "\n",
    "kappa_list = [1e-6, 1, 1e6]\n",
    "### End of parameter definitions\n",
    "\n",
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "mbda = 1\n",
    "\n",
    "def chi(x, y):\n",
    "    return np.logical_and(np.greater(x, 0.5), np.greater(y, 0.5))\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "ux = sympy.sin(2*pi * x) * sympy.sin(2 * pi * y)\n",
    "uy = sympy.cos(pi * x) * (y-0.5)**2\n",
    "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "divu = dux_x + duy_y\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n",
    "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n",
    "\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = np.zeros((num_refs, len(kappa_list)))\n",
    "    flux_err = np.copy(u_err)\n",
    "    \n",
    "    for iter1, g in enumerate(setup_grids.grid_sequence_fixed_lines(basedim, num_refs, grid_type, pert, subdom_func=chi)):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        for iter2, kappa in enumerate(kappa_list):        \n",
    "        # Permeability tensor\n",
    "            char_func_cells = chi(g.cell_centers[0], g.cell_centers[1]) * 1.\n",
    "            mat_vec = (1 - char_func_cells) + kappa * char_func_cells\n",
    "\n",
    "            k = stiffness(2, mat_vec, mat_vec)\n",
    "\n",
    "            # Set type of boundary conditions - Dirichlet\n",
    "            bound_faces = g.get_all_boundary_faces()\n",
    "            bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "\n",
    "            # MPFA discretization, and system matrix\n",
    "            stress, bound_stress = mpsa.mpsa(g, k, bound_cond)\n",
    "            div = fvutils.vector_divergence(g)\n",
    "            a = div * stress\n",
    "\n",
    "            # Boundary conditions\n",
    "            xf = g.face_centers\n",
    "            char_func_bound = chi(xf[0, bound_faces], xf[1, bound_faces]) * 1\n",
    "            u_bound = np.zeros((g.dim, g.num_faces))\n",
    "            u_bound[0, bound_faces] = ux_f(xf[0, bound_faces], xf[1, bound_faces])\\\n",
    "                    / ((1 - char_func_bound) + kappa * char_func_bound)\n",
    "            u_bound[1, bound_faces] = uy_f(xf[0, bound_faces], xf[1, bound_faces]) \\\n",
    "                    / ((1 - char_func_bound) + kappa * char_func_bound)\n",
    "\n",
    "            # Right hand side - contribution from the solution and the boundary conditions\n",
    "            xc = g.cell_centers\n",
    "            rhs = np.vstack((rhs_x_f(xc[0], xc[1]), \n",
    "                             rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n",
    "            b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "            # Solve system, derive fluxes\n",
    "            u_num = spsolve(a, b)\n",
    "            stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "\n",
    "            ux_num = u_num[::2]\n",
    "            uy_num = u_num[1::2]\n",
    "            \n",
    "            stress_x_num = stress_num[::2]\n",
    "            stress_y_num = stress_num[1::2]\n",
    "            # Exact solution\n",
    "            ux_ex = ux_f(xc[0], xc[1])/ ((1 - char_func_cells) + kappa * char_func_cells)\n",
    "            uy_ex = uy_f(xc[0], xc[1])/ ((1 - char_func_cells) + kappa * char_func_cells)\n",
    "            u_ex = np.vstack((ux_ex, uy_ex))\n",
    "            u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "\n",
    "            sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n",
    "            sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n",
    "\n",
    "            stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n",
    "            stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n",
    "            stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                     stress_y_num - stress_y_ex))\n",
    "            stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n",
    "            \n",
    "            u_err[iter1, iter2] = np.sqrt(np.sum(g.cell_volumes * u_diff**2)) \\\n",
    "                        / np.sqrt(np.sum(g.cell_volumes * u_ex**2))\n",
    "            flux_err[iter1, iter2] = np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))\\\n",
    "                        / np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2))\n",
    "\n",
    "    return u_err, flux_err\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(f_cart_nopert,'stress error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'displacement error')\n",
    "print(f_cart_pert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(f_triang_nopert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('triangular', 0.5)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'displacement error')\n",
    "print(f_triang_pert,'stress error')\n",
    "\n",
    "# EK: These values were hard-coded 12.06.2017.\n",
    "assert np.abs(u_cart_nopert[2, 0] - 0.01980355  ) < 1e-7\n",
    "assert np.abs(f_cart_nopert[2, 0] - 0.02181871   ) < 1e-7\n",
    "assert np.abs(u_cart_pert[2, 0] - 0.02294892  ) < 1e-7\n",
    "assert np.abs(f_cart_pert[2, 0] - 0.02941291  ) < 1e-7\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2, 0] - 0.01574973  ) < 1e-7\n",
    "assert np.abs(f_triang_nopert[2, 0] - 0.01890709  ) < 1e-7\n",
    "assert np.abs(u_triang_pert[2, 0] - 0.01911508  ) < 1e-7\n",
    "assert np.abs(f_triang_pert[2, 0] - 0.02295011   ) < 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Incompressible limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[[0.03697053 0.03753405 0.03759228]\n",
      " [0.00976093 0.00987151 0.00988294]\n",
      " [0.00248737 0.00251144 0.00251393]\n",
      " [0.00062543 0.00063108 0.00063166]] displacement error\n",
      "[[0.05751235 0.05696182 0.05690496]\n",
      " [0.0198125  0.01970924 0.01969853]\n",
      " [0.00646923 0.00645496 0.00645348]\n",
      " [0.00209802 0.00209683 0.0020967 ]] stress error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[[0.03956063 0.04032235 0.04040406]\n",
      " [0.01064916 0.01079559 0.01081118]\n",
      " [0.00270644 0.00277624 0.00279648]\n",
      " [0.0007012  0.00071505 0.00072281]] displacement error\n",
      "[[0.10770594 0.54796429 5.20910705]\n",
      " [0.03973921 0.24607287 2.44840796]\n",
      " [0.01322056 0.02545624 0.19096982]\n",
      " [0.00524277 0.0151903  0.18131059]] stress error\n",
      "Triangular errors - no perturbations\n",
      "[[0.01757125 0.01776094 0.01778063]\n",
      " [0.00429797 0.00432499 0.00432781]\n",
      " [0.00108379 0.00108757 0.00108797]\n",
      " [0.00027397 0.00027454 0.0002746 ]] displacement error\n",
      "[[0.033668   0.03356879 0.03355879]\n",
      " [0.01149681 0.0115022  0.01150278]\n",
      " [0.00390429 0.00391628 0.00391752]\n",
      " [0.00132659 0.00133252 0.00133313]] stress error\n",
      "Triangular errors - perturbation 0.5\n",
      "[[0.02128916 0.0216464  0.0216837 ]\n",
      " [0.00499226 0.00503323 0.00503754]\n",
      " [0.00123134 0.00123924 0.00124007]\n",
      " [0.00030901 0.00031037 0.00031051]] displacement error\n",
      "[[2.15700712e-01 2.14964132e+00 2.15497574e+01]\n",
      " [2.19628718e-02 8.79909686e-02 8.55960998e-01]\n",
      " [7.91132377e-03 1.15770525e-02 8.16011471e-02]\n",
      " [3.26020037e-03 5.74393379e-03 4.76682084e-02]] stress error\n"
     ]
    }
   ],
   "source": [
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "lmbda = 1\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 4\n",
    "grid_type = 'cart'\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "ux = sympy.sin(y) * (sympy.cos(x)-1)\n",
    "uy = sympy.sin(x) * (1-sympy.cos(y))\n",
    "ux_f = sympy.lambdify((x, y), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y), uy, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "divu = dux_x + duy_y\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y), sxy, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y), syy, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y)\n",
    "rhs_x_f = sympy.lambdify((x, y), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y), rhs_y, 'numpy')\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = np.zeros((num_refs, 3))\n",
    "    flux_err = np.zeros((num_refs, 3))\n",
    "\n",
    "    for iter1, g in enumerate(setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain)):\n",
    "        for iter2, exponent in enumerate(range(2, 5)):\n",
    "            # Reset the random seed for every grid realization.\n",
    "            # This should make no difference for the convergence test, \n",
    "            # but it makes sure that we can run unit tests based on the values obtained\n",
    "            # here.\n",
    "            np.random.seed(42)\n",
    "\n",
    "            lmbda = 10**exponent\n",
    "\n",
    "            # Permeability tensor\n",
    "            mu_c = mu * np.ones(g.num_cells)\n",
    "            lmbda_c = lmbda * np.ones(g.num_cells)\n",
    "            k = stiffness(2, mu_c, lmbda_c)\n",
    "\n",
    "            # Set type of boundary conditions - Dirichlet\n",
    "            bound_faces = g.get_all_boundary_faces()\n",
    "            bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "\n",
    "            # MPFA discretization, and system matrix\n",
    "            stress, bound_stress = mpsa.mpsa(g, k, bound_cond)\n",
    "            div = fvutils.vector_divergence(g)\n",
    "            a = div * stress\n",
    "\n",
    "            # Boundary conditions\n",
    "            xf = g.face_centers\n",
    "            u_bound = np.zeros((g.dim, g.num_faces))\n",
    "            u_bound[0, bound_faces] = ux_f(xf[0, bound_faces], xf[1, bound_faces])\n",
    "            u_bound[1, bound_faces] = uy_f(xf[0, bound_faces], xf[1, bound_faces])\n",
    "\n",
    "\n",
    "            # Right hand side - contribution from the solution and the boundary conditions\n",
    "            xc = g.cell_centers\n",
    "            rhs = np.vstack((rhs_x_f(xc[0], xc[1]), rhs_y_f(xc[0], xc[1]))) * g.cell_volumes\n",
    "            b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "            # Solve system, derive fluxes\n",
    "            u_num = spsolve(a, b)\n",
    "            stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "\n",
    "            ux_num = u_num[::2]\n",
    "            uy_num = u_num[1::2]\n",
    "\n",
    "            stress_x_num = stress_num[::2]\n",
    "            stress_y_num = stress_num[1::2]\n",
    "            # Exact solution\n",
    "            ux_ex = ux_f(xc[0], xc[1])\n",
    "            uy_ex = uy_f(xc[0], xc[1])\n",
    "            u_ex = np.vstack((ux_ex, uy_ex))\n",
    "            u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "\n",
    "            sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1]), sxy_f(xf[0], xf[1])))\n",
    "            sy_ex_faces = np.vstack((syx_f(xf[0], xf[1]), syy_f(xf[0], xf[1])))\n",
    "\n",
    "            stress_x_ex = np.sum(g.face_normals[:2] * sx_ex_faces, axis=0)\n",
    "            stress_y_ex = np.sum(g.face_normals[:2] * sy_ex_faces, axis=0)\n",
    "            stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                     stress_y_num - stress_y_ex))\n",
    "            stress_ex = np.vstack((stress_x_ex, stress_y_ex))\n",
    "        \n",
    "            u_err[iter1, iter2] = np.sqrt(np.sum(g.cell_volumes * u_diff**2)) \\\n",
    "                        / np.sqrt(np.sum(g.cell_volumes * u_ex**2))\n",
    "            flux_err[iter1, iter2] = np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))\\\n",
    "                        / np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2))\n",
    "\n",
    "    return u_err, flux_err\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(f_cart_nopert,'stress error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'displacement error')\n",
    "print(f_cart_pert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(f_triang_nopert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('triangular', 0.5)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'displacement error')\n",
    "print(f_triang_pert,'stress error')\n",
    "\n",
    "# EK: These values were hard-bcoded 12.06.2017.\n",
    "assert np.abs(u_cart_nopert[2, 0] -  0.00248737   ) < 1e-7\n",
    "assert np.abs(f_cart_nopert[2, 0] -    0.00646923  ) < 1e-7\n",
    "assert np.abs(u_cart_pert[2, 0] -  0.00270644   ) < 1e-7\n",
    "assert np.abs(f_cart_pert[2, 0] -  0.01322056   ) < 1e-7\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2, 0] -  0.00108379   ) < 1e-7\n",
    "assert np.abs(f_triang_nopert[2, 0] -   0.00390429  ) < 1e-7\n",
    "assert np.abs(u_triang_pert[2, 0] -   0.00123134   ) < 1e-7\n",
    "assert np.abs(f_triang_pert[2, 0] -   7.91132377e-03 ) < 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that on triangular grids, the convergence order for stress deteriorates for high Poisson's ratios. See Keilegavlen, Nordbotten IJNME 2017 for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3d homogeneous problem\n",
    "\n",
    "This example also compares computational times for discretization and two linear solvers (direct and amg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time for discretization 1.08636474609375\n",
      "Time for AMG 0.020662784576416016\n",
      "Time for discretization 1.6143503189086914\n",
      "Time for AMG 0.07893586158752441\n",
      "Time for discretization 1.8830769062042236\n",
      "Time for AMG 0.4144134521484375\n",
      "Cartesian errors - no perturbations\n",
      "[0.14582494480280667, 0.036685508143867125, 0.009623477148373078] displacement error\n",
      "[0.14173578960074407, 0.05453461058511372, 0.018777567804334296] stress error\n",
      "Time for discretization 1.1395857334136963\n",
      "Time for AMG 0.010153055191040039\n",
      "Time for discretization 1.394442081451416\n",
      "Time for AMG 0.0783224105834961\n",
      "Time for discretization 1.727994680404663\n",
      "Time for AMG 0.37769436836242676\n",
      "Cartesian errors - perturbation 0.5\n",
      "[0.1823425575505196, 0.03669658584037131, 0.009785391062793553] displacement error\n",
      "[0.15394769153699012, 0.06655533666156926, 0.024831033744125167] stress error\n",
      "Time for discretization 1.5092995166778564\n",
      "Time for AMG 0.04783058166503906\n",
      "Time for discretization 2.5991227626800537\n",
      "Time for AMG 0.35640907287597656\n",
      "Time for discretization 5.721246719360352\n",
      "Time for AMG 4.348304986953735\n",
      "Triangular errors - no perturbations\n",
      "[0.06256710001786384, 0.015518298189089326, 0.004112646550453208] displacement error\n",
      "[0.09819533736107755, 0.040432049644100486, 0.01608810670167512] stress error\n",
      "Time for discretization 1.1475048065185547\n",
      "Time for AMG 0.04819488525390625\n",
      "Time for discretization 1.789534568786621\n",
      "Time for AMG 0.342240571975708\n",
      "Time for discretization 6.363281965255737\n",
      "Time for AMG 4.3170764446258545\n",
      "Triangular errors - perturbation 0.5\n",
      "[0.07257865039340129, 0.01699696434861487, 0.004336430438721095] displacement error\n",
      "[0.10603237914286788, 0.04292864828702438, 0.017041050867287796] stress error\n"
     ]
    }
   ],
   "source": [
    "# Permeability tensor, scalar for simplicity\n",
    "mu = 1\n",
    "lmbda = 1\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 2\n",
    "domain = np.array([1, 1, 1])\n",
    "basedim = np.array([base, base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "# Analytical solution. Essentially randomly made\n",
    "x, y, z = sympy.symbols('x y z')\n",
    "ux = sympy.sin(x) * sympy.cos(y) * (0.5 - z)**2\n",
    "uy = sympy.sin(y) * x**2 * sympy.cos(z)\n",
    "uz = sympy.cos(x) * sympy.sin(z) * y * (1-z)\n",
    "ux_f = sympy.lambdify((x, y, z), ux, 'numpy')\n",
    "uy_f = sympy.lambdify((x, y, z), uy, 'numpy')\n",
    "uz_f = sympy.lambdify((x, y, z), uz, 'numpy')\n",
    "dux_x = sympy.diff(ux, x)\n",
    "dux_y = sympy.diff(ux, y)\n",
    "dux_z = sympy.diff(ux, z)\n",
    "duy_x = sympy.diff(uy, x)\n",
    "duy_y = sympy.diff(uy, y)\n",
    "duy_z = sympy.diff(uy, z)\n",
    "duz_x = sympy.diff(uz, x)\n",
    "duz_y = sympy.diff(uz, y)\n",
    "duz_z = sympy.diff(uz, z)\n",
    "\n",
    "divu = dux_x + duy_y + duz_z\n",
    "\n",
    "sxx = 2 * mu * dux_x + lmbda * divu\n",
    "sxy = mu * (dux_y + duy_x)\n",
    "sxz = mu * (dux_z + duz_x) \n",
    "syx = mu * (duy_x + dux_y)\n",
    "syy = 2 * mu * duy_y + lmbda * divu\n",
    "syz = mu * (duy_z + duz_y)\n",
    "szx = mu * (duz_x + dux_z)\n",
    "szy = mu * (duy_z + duz_y)\n",
    "szz = 2 * mu * duz_z + lmbda * divu\n",
    "\n",
    "sxx_f = sympy.lambdify((x, y, z), sxx, 'numpy')\n",
    "sxy_f = sympy.lambdify((x, y, z), sxy, 'numpy')\n",
    "sxz_f = sympy.lambdify((x, y, z), sxz, 'numpy')\n",
    "syx_f = sympy.lambdify((x, y, z), syx, 'numpy')\n",
    "syy_f = sympy.lambdify((x, y, z), syy, 'numpy')\n",
    "syz_f = sympy.lambdify((x, y, z), syz, 'numpy')\n",
    "szx_f = sympy.lambdify((x, y, z), szx, 'numpy')\n",
    "szy_f = sympy.lambdify((x, y, z), szy, 'numpy')\n",
    "szz_f = sympy.lambdify((x, y, z), szz, 'numpy')\n",
    "\n",
    "rhs_x = sympy.diff(sxx, x) + sympy.diff(syx, y) + sympy.diff(szx, z)\n",
    "rhs_y = sympy.diff(sxy, x) + sympy.diff(syy, y) + sympy.diff(szy, z)\n",
    "rhs_z = sympy.diff(sxz, x) + sympy.diff(syz, y) + sympy.diff(szz, z)\n",
    "rhs_x_f = sympy.lambdify((x, y, z), rhs_x, 'numpy')\n",
    "rhs_y_f = sympy.lambdify((x, y, z), rhs_y, 'numpy')\n",
    "rhs_z_f = sympy.lambdify((x, y, z), rhs_z, 'numpy')\n",
    "\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "            \n",
    "        # Permeability tensor\n",
    "        mu_c = mu * np.ones(g.num_cells)\n",
    "        lmbda_c = lmbda * np.ones(g.num_cells)\n",
    "        k = stiffness(g.dim, mu_c, lmbda_c)\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.get_boundary_faces()\n",
    "        bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        tm = time()\n",
    "        stress, bound_stress = mpsa.mpsa(g, k, bound_cond)\n",
    "        print('Time for discretization ' + str(time() - tm))\n",
    "        div = fvutils.vector_divergence(g)\n",
    "        a = div * stress\n",
    "        a.tocsr()\n",
    "        \n",
    "        # Boundary conditions\n",
    "        xf = g.face_centers\n",
    "        u_bound = np.zeros((g.dim, g.num_faces))\n",
    "        u_bound[0, bound_faces] = ux_f(xf[0, bound_faces], xf[1, bound_faces], xf[2, bound_faces])\n",
    "        u_bound[1, bound_faces] = uy_f(xf[0, bound_faces], xf[1, bound_faces], xf[2, bound_faces])\n",
    "        u_bound[2, bound_faces] = uz_f(xf[0, bound_faces], xf[1, bound_faces], xf[2, bound_faces])\n",
    "        \n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        rhs = np.vstack((rhs_x_f(xc[0], xc[1], xc[2]), \n",
    "                         rhs_y_f(xc[0], xc[1], xc[2]),\n",
    "                         rhs_z_f(xc[0], xc[1], xc[2]))) * g.cell_volumes\n",
    "        b = rhs.ravel('F') - div * bound_stress * u_bound.ravel('F')\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        tm = time()\n",
    "        smooth = smooth=('energy', {'krylov': 'gmres', 'degree': 2})\n",
    "        \n",
    "        null_space = np.zeros((g.dim * g.num_cells, 6))\n",
    "        null_space[0::g.dim, 0] = 1\n",
    "        null_space[1::g.dim, 1] = 1\n",
    "        null_space[2::g.dim, 2] = 1\n",
    "        # Rotations\n",
    "        null_space[0::g.dim, 3] = 1\n",
    "        null_space[1::g.dim, 3] = -1\n",
    "        null_space[1::g.dim, 4] = 1\n",
    "        null_space[2::g.dim, 4] = -1\n",
    "        null_space[2::g.dim, 5] = 1\n",
    "        null_space[0::g.dim, 5] = -1\n",
    "\n",
    "        res = []\n",
    "        a = a.tocsr()\n",
    "        ml = pyamg.smoothed_aggregation_solver(a, null_space, smooth=smooth, symmetry='nonsymmetric')\n",
    "        u_num = ml.solve(b, residuals=res, tol=1e-12)\n",
    "        stress_num = stress * u_num + bound_stress * u_bound.ravel('F')\n",
    "        print('Time for AMG ' + str(time() - tm))\n",
    "        ux_num = u_num[::g.dim]\n",
    "        uy_num = u_num[1::g.dim]\n",
    "        uz_num = u_num[2::g.dim]\n",
    "        \n",
    "        stress_x_num = stress_num[::g.dim]\n",
    "        stress_y_num = stress_num[1::g.dim]\n",
    "        stress_z_num = stress_num[2::g.dim]\n",
    "        # Exact solution\n",
    "        ux_ex = ux_f(xc[0], xc[1], xc[2])\n",
    "        uy_ex = uy_f(xc[0], xc[1], xc[2])\n",
    "        uz_ex = uz_f(xc[0], xc[1], xc[2])\n",
    "        u_ex = np.vstack((ux_ex, uy_ex, uz_ex))\n",
    "        u_diff = np.vstack((ux_num - ux_ex, uy_num - uy_ex))\n",
    "        \n",
    "        sx_ex_faces = np.vstack((sxx_f(xf[0], xf[1], xf[2]),\n",
    "                                 sxy_f(xf[0], xf[1], xf[2]),\n",
    "                                 sxz_f(xf[0], xf[1], xf[2])))\n",
    "        \n",
    "        sy_ex_faces = np.vstack((syx_f(xf[0], xf[1], xf[2]),\n",
    "                                 syy_f(xf[0], xf[1], xf[2]),\n",
    "                                 syz_f(xf[0], xf[1], xf[2])))\n",
    "        \n",
    "        sz_ex_faces = np.vstack((szx_f(xf[0], xf[1], xf[2]),\n",
    "                                 szy_f(xf[0], xf[1], xf[2]),\n",
    "                                 szz_f(xf[0], xf[1], xf[2])))\n",
    "\n",
    "        stress_x_ex = np.sum(g.face_normals[:g.dim] * sx_ex_faces, axis=0)\n",
    "        stress_y_ex = np.sum(g.face_normals[:g.dim] * sy_ex_faces, axis=0)\n",
    "        stress_z_ex = np.sum(g.face_normals[:g.dim] * sz_ex_faces, axis=0)\n",
    "        stress_diff = np.vstack((stress_x_num - stress_x_ex, \n",
    "                                 stress_y_num - stress_y_ex,\n",
    "                                 stress_z_num - stress_z_ex))\n",
    "        stress_ex = np.vstack((stress_x_ex, stress_y_ex, stress_z_ex))\n",
    "        \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * u_diff**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * stress_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * stress_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'displacement error')\n",
    "print(f_cart_nopert,'stress error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'displacement error')\n",
    "print(f_cart_pert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('tetrahedral', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'displacement error')\n",
    "print(f_triang_nopert,'stress error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('tetrahedral', 0.3)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'displacement error')\n",
    "print(f_triang_pert,'stress error')\n",
    "\n",
    "# EK: These values were hard-coded 16.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] - 0.0096234771484497109) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.018777567804377917) < 1e-10\n",
    "assert np.abs(u_cart_pert[2] - 0.0097853910627935265) < 1e-10\n",
    "assert np.abs(f_cart_pert[2] - 0.024831033744125156) < 1e-10\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2] - 0.004112646550453208) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.01608810670167512) < 1e-10\n",
    "assert np.abs(u_triang_pert[2] - 0.004336430438721095) < 1e-10\n",
    "assert np.abs(f_triang_pert[2] - 0.017041050867287796) < 1e-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyamg.smoothed_aggregation_solver?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
